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Abstract 


This study presents an algorithm for solving optimal control problems with 
the objective function of the Lagrange-type and multiple delays on both the 
state and control variables of the constraints, with bounds on the control 
variable. The full discretization of the objective functional and the multi- 
ple delay constraints is carried out by using the Simpson numerical scheme. 
The discrete recurrence relations generated from the discretization of both 
the objective functional and constraints are used to develop the matrix op- 
erators, which satisfy the basic spectral properties. The primal-dual residu- 
als of the algorithm are derived in order to ascertain the rate of convergence 
of the algorithm, which performs faster when relaxed with an accelerator 
variant in the sense of Nesterov. The direct numerical approach for han- 
dling the multi-delay control problem is observed to obtain an accurate 
result at a faster rate of convergence when over-relaxed with an accelerator 
variant. This research problem is limited to linear constraints and objec- 
tive functional of the Lagrange-type and can address real-life models with 
multiple delays as applicable to quadratic optimization of intensity mod- 
ulated radiation theory planning. The novelty of this research paper lies 
in the method of discretization and its adaptation to handle linearly and 
proximal bound-constrained program formulated from the multiple delay 
optimal control problems. 
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1 Introduction 


Optimal control problems (OCPs) with delays in state and/or control vari- 
ables play important roles in the modeling of real-life scenarios. Time delay 
effects in economical, physical, chemical, and biological processes involving 
optimal control systems cannot be neglected, especially when they involve 
the transmission of information between different parts of the system. These 
hereditary effects of constant delays have been deliberated upon by Xue and 
Duanin [27], Rihan and Anwar [21], Laarabi, Abta, and Hattaf, [14] and 
Bashier and Patidar [1]. Gollman, Kern and Maurer [9] did extensive work 
on mixed control-state inequality constraints with single delay on both state 
and delay variables with the initial and terminal boundary conditions in a 
general mixed form. Later, Gollmann and Maurer [10] extended their work 
to multiple time delays in the control and state variables with mixed control- 
state constraints. Olotu and Dawodu [18] and Olotu, Dawodu, and Yusuf [19] 
worked on one-dimensional and generalized control problems, respectively, 
with multiple delay constants and unbounded control variable (vector) using 
modified alternating direction method of multipliers (ADMM), while their 
earlier work in [17] on the delay proportional control system was done by the 
method of Quasi-Newton embedded augmented Lagrangian functional. 

The Douglas—Rachford (D-R) splitting method called the ADMM (AI- 
ternating Direction Method of Multipliers) had been deployed by a plethora 
of authors in solving OCPs; though very few have deployed it to solving 
multi-delay OCPs with control bounds. The patronage is based on the fact 
that ADMM enjoys the strong convergence properties of the method of mul- 
tipliers and the decomposability property of dual ascent as in the works of 
Boley [2], Boyd and Vandenberghe [4], Sun, Toh, and Yang [23], and Yao et. 
al [28]; hence its relevance in solving decentralized and block-convex opti- 
mization problems like the large-scaled power or distributed systems cannot 
be overemphasized. Yao et al. [28] discovered that the ADMM could be 
used as a tool in solving a vehicle routing problem structured as an inte- 
ger programming problem provided; the quadratic penalty terms used in 
ADMM can be reduced into simple linear functions. He and Yuan [13] and 
He, Hong, and Yuan [12] worked on the proximal Jacobian decomposition 
of the augmented Lagrangian method for the multiple-block separable con- 
vex minimization problem. It was observed that it has a strong affinity for 
ADMM. 

Many authors have researched constrained convex quadratic problems 
with simple bounds. Voglis and Lagaris [24] was among the early authors 
that developed an algorithm for quadratic programming problems with box 
constrained using the active set strategy with applications mainly to circus 
tent, random and bi-harmonic problems. However, the scalability, simplicity, 
and potential of the ADMM in solving large-scaled structured convex opti- 
mization problems made it widely acceptable. The ADMM has the ability to 
solve the primal-dual feasibility updates in parallel and is well amenable to 
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the Gauss-Seidel accelerator or relaxation in order to improve its efficiency 
and rate of convergence as in Nesterov [16]. 

Pasquale and Gerardo [20] also reviewed the main results on global op- 
timality conditions for establishing the global minimization of a quadratic 
problem with box constraints. However, Carreira-Perpindn [5] worked on 
the proximal bound-constrained quadratic problem subject to simple box 
constraints. The ADMM algorithm was considered using the convex proxi- 
mal operator by Moreau [15], Rockafellar [22], and Combettes and Pesquet 
in [6]. The algorithm was found to be particularly efficient in solving a collec- 
tion of proximal operators that share a same quadratic form with a high rate 
of convergence and accuracy. Wu and Shang [25] further presented a global 
optimization method for solving general nonlinear programming problems 
subject to box constraints using a differential flow on the dual feasible space. 
The criteria for the global optimality and existence of complete solutions to 
the original problems were also given. 

Xinmin et al. [26] researched the application of the ADMM algorithm 
to the constrained quadratic optimization of intensity modulated radiation 
theory planning. A graph form convex optimization model was formulated 
as minh(y) such that y= Az, 0 <a < ug, ly < x < uy, where the meaning 
and conditions for the bounds were well-defined. This approach involves the 
construction of a proximal operator and the formulation of the optimization 
problem in the graph form of the ADMM algorithm. In the implementation 
of the algorithm, a clinically acceptable solution was obtained with fewer 
memory requirements when compared with other optimization solvers. 


2 Statement of problem 


The generalized form of the multiple delay OCP with bounded (box) control is 
expressed below with quadratic continuous functional F and linear dynamical 
function f. 


2.1 Optimal control model 


Considering the multi-delay OCP with bounded control below: 


T 
min J(x, u) =; / F(t, x(t), u(t)) dt (1) 
d é 
s.t. &(t) <Ax + But+ 5° aja(t— r;) +S° Bult — qi) (2) 
j= l=1 
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s.t. &(t) <f(t,2(t),u(t),c*(t—r;),u"(t—q)), 5) =1,2,...,d(e), (3) 

x(to) =Xo, (4) 

a(t)=9(t), to—r<t<to, (5) 

u(t) =p(t),  to—q<t<to, (6) 

¥ < ult) <o, (7) 
where 


a Te 
— Daiji anr eee Oy) € R", = (61, b2,---,n) € R", 
U1, U2,--- ia) € R”, 


( 
( 
z"(t —7rj) = (a(t —11), a(t —re),...,2(t—ra)) € R®, 
pa 
( 
( 


u(t — qi), u(t — q2),..-,u(t—qe)) ER™, 
1, W2,...,0m)? €R™, 

V1y.-65Y%m)? € R™, 

= (01,...,0m)’ € R™, 


r= max{rj}4_1, and g = max{q}/_1. 


2.2 Assumptions and properties of the model 


The underlining assumptions and properties of the functions are stated below: 
(i) The time interval I = [to, T] C R is a subset of real numbers R; 


(ii) The state function x(t) : I > X of the optimal control system is abso- 
lutely continuous on I provided the subset X is closed and bounded in 
R”: 

(iii) The control function u(t) : I > U of the control system is piece-wise 


continuous and Lebesgue measurable on I provided the subset U is 
closed and bounded in R™; 


(iv) The numbers rj,7r2,...,rq and qi, q2,---,;Ge are the known monotoni- 
cally increasing positive delay constants for the state and control vari- 
ables, respectively; that is, rj <1j41 for all r; € R for 7 = 1,2,...,d 
and q < qi41 for all q € R for 1 =1,2,...,¢; 


(v) The initial (pre-shaped) functions p(t) : [to — r, to] 4 R” and w(t) : 
[to — g,to] + R™ for the state and control variables, respectively, are 
known and piecewise continuous; 


(vi) The numbers r = max{r; en and q = max{q}/_, represent the values 
of the maximum delay constants for the state and control variables, 
respectively; 
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(vii) The cost functional F' : [to,T] x R” x R™ — RU {oo} is nonlinear, 
piecewise continuous and measurable; 


(viii) The constraint functional f : [to, T]x RO+9” x RU+9™ _, R” is linear 
and piece-wise continuous; 


(ix) There exist an admissible pair p = (a(-),u(-)) € Q := ([to, T],R”) x 
([to, T],R”) that satisfies (1)-(7) of the multi-delay OCP, where the 
functional J(x(-),u(-)) is minimized by the the admissible triple Q := 


(x(-), u-), PC). 


3 Methodology 


A very simple discretization of (1) and (3) is obtained by means of the Simp- 
sons and Euler numerical schemes, respectively. Recurrence relations gen- 
erated are used in developing the respective large sparse matrix operators. 
The augmented Lagrangian functional with relaxation parameter (factor) is 
then used as an accelerator variant in the formulation of the modified ADMM 
(M-ADMM) algorithm. The spectral and convergence analyses are carried 
out to ensure the well-posedness of the algorithm. 


3.1 Background and preliminaries 


Suppose that the discrete time interval is given as I, = [tx, te4i] by letting 
t, = to + kd for k = 0,1,...,N. Then the discretization operator, say fz, 
maps each discrete point t, in the discrete time interval J, C R into each 
discrete (grid) point of the concatenated state vector a) € R™ for all i= 
1,2,...,n, while the operator, say f,,, maps the points into each discrete point 
of the concatenated control vector ul’) € R™A4+) for all j =1,2,...,m. It 
is then expressed as 


fu : [to, T] Cc R— iu € U5 (9) 
where = (&1,...,2n) € R™ and @ = (t,...,i%m) € R™ +t” for 
ti = (a\?,..., 06%) € RN, and a; = GP eos) € RN+1, respec- 


tively. Considering the existence of multiple delay constants r; € R (for 
j =1,2,...,d) and q € R (for / = 1,2,...,e) on the state and control vari- 
ables, respectively. Then some basic Theorems are expressed below to aid in 
the process of discretization. 
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Theorem 1 (Rationality theorem). Given the real numbers rj, 7341 > 0 for 
rj <1j41, then there exists a unique real number 6 < 1 such that the ratios 
of the numbers is a rational number Q. 


Proof. Let the delay constants be arranged in an increasing manner; that is, 
rj <7rj41 and qm < qm4+1 for every j, 1. It then implies that for every r; and 
qm, there exist integers m; € Z* and q € Zt such that 
} 

ge Ege (10) 

™4 Ow] 
where 6 is known as the step-length or interval length. This clearly shows 
that the delays, for both the state and control variables, are integer multiples 
of the step-length. 


The proof of this transformation technique was suggested by Guinn [11] 
to derive first-order necessary conditions for unconstrained OCPs with pure 
state delays. To further elaborate on the assumptions above and the appli- 
cations, we then introduce Theorem 2 below. 


Theorem 2. Given any interval [a,b], there exists a step-length h € Rt such 
that each sub-interval is a constant multiple of the step-length. 


Proof. Let the entire length of the interval 1 = b — a be partitioned into m 
sub-intervals [rj1,7r;] for 7 = 1,2,...,m. Then there exists h € R with 
r; = hv; for all v; € Z* such that the length of each sub-interval can be 
written as L; = 1; —Tj-1- 


by _ (ry = Tj-1) 


pa SPE a nj + 6; nj €Z*, 6; € (0,1). 
Then, 
hu; — ho;_ 
( 7 h 7 1) = (vu; 05-4) => Nh; +6; 5 E Z*, 5; = [0, 1); (11) 


Since v; — vj;-1 € Z* then n; + 6; € Zt such that 6; = 0 or 6; = 1 for all 
j =1,2,...,d (by Theorem 1). Since 6; = 1 ¢ [0,1), then 6; =0 for all 7 = 
1,2,...,d. By slotting 6; = 0 into (11), we then conclude that the length 
of each sub-interval 1; = (rj; — rj-1) can be expressed as a multiple integer 
vj € Z* of h expressed as rj = 7;7-1 + hn;, where vj; — vj-1 = nj for all 
9 Hg Doe gM: 


The essence of Theorems 1 and 2 above is to establish the fact that there 
exists a real number, referred to as step-length, h € R that divides the entire 
state and control delay constants with multiples of real positive integers. 


Consequent to 1 and 2, the multiple state and control variables are pre- 
sented as follows: 
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P(tk—v;) k—vj; <0; k=0,1,2,...,v;, t € [to —r,to], (known) 
ak) =o glk), kv; > 0; k= (vj +1),.-.,N, t € [to,T], (unknown) (12) 
x(to) = zo given 


and 


or ees k—w, <0; k=0,1,2,...,(w;—1), t € [to —4, to], (known) 
ulk-“), «kk —w, > 0 for k = wy, (w; +1),...,N, t € [to, T], (unknown) 
(13) 
where r = rq = max{rj}4_4 and q = ge = max{q}j_, are the largest 
delays on the state and control variables, respectively. In other words, the 
discretized delay state and control variables are represented in the form below: 


P 
ata) (0%) a9), ao)" ERY, forany (ld) 
T 
yh) = Cig oe a) ER”, for any J, (15) 


while x(t) and u(t) are estimated within their respective delay intervals by 
the given delay functions ¢(t) and w(t), respectively, such that (t) ~ (t_x) 
and w(t) ~ W(t_x), where t, = to — kd for k = 1,2,...,va,..., We if va < We. 
However, the largest state and control delay vectors for each iterate k = 
0,1,... are represented below: 


T 
*= (x0), 2 (2-03) ga) ER™Y-*) for all j = 1,2,...,, 


(16) 
T 
a= (ui, aig Same a) | Ee R™%N+1) for all | =1,2,...,. 
(17) 
3.2 Discretization of the objective functional 
The objective function in (1) is of the form 
1 T 
min J(z,u) = >| (a? Px + u™ Qu) dt, (18) 
to 
where z = (21,22,.-.,2n)’ € R”, u = (a1, ta,...,Um)? € R™, P € R"™”, 


and Q@ € R™*™. The full discretization of the objective function (1) using 
the Simpson integration formula gives 


i 
7 y PG) 
min J (2, u) SS F(tg, 2 ) (19) 
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and when expanded, it can be expressed as the quadratic function 


1 = 1 = 
min 52° Pr + 50" Qu+t R, (20) 
where R = 2 (¢(0))T Pe € R and the concatenated state and control vari- 
ables are & = (w) ,2,...,a)? eR" and a = (u,u,...,uQ)T € 


RN 7a) respectively. However, the block-diagonal coefficient matrices 
P € R°N*N and Q € R™N+DxmWN+) are the block-matrix operators 
of the objective functional stated below as 


5 
GO i). 2% 414-282 
48P Q wee vas 0 e 
45 
Q 2p De 
7 3 7 a 2G” 
LS and Q= 3 
SP 0 
3 
ee 0 §P 2Q 0 
(QV) ale ale ade 0 38 


The detailed derivation is in Appendix as in Olotu, Dawodu, and Yusuf [19]. 
However, in determining the properties of the derived matrix operators above, 
the following definitions are then introduced. 


Definition 1 (Sylvester criterion). A square (nxn) matrix is positive definite 
if it is real, symmetric with all the eigenvalues (A; > 0; i = 1,2,...,n) or the 
leading principal minors (M; > 0; i =1,2,...,n) are strictly positive. 


Definition 2 (Corollary to Sylvester criterion). If the principal diagonal 
entries a;; are the only nonzero entries of a square (n x n) matrix such that 
a;; = 0 for i ~ j, then the leading principal minors (M; > 0; i = 1,2,...,n) 
are given by M; = me ag tort = 1,2.4045% 


Definition 3. Every positive definite matrix is invertible. 


Sequel to Definitions 1 and 2, the constructed matrix-operator P is symmetric 
and positive definite since all the leading principal minors are strictly positive. 
Consequently, by Definition 3, the matrix P is invertible. 


3.3 Discretization of the constraints 


The Euler numerical scheme in (21) below is used for the discretization of 
the multi-delay constraints; 


ohh) — gl®) 4 gf, (21) 
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where f* = f(t, 2, u”,#™, a) and expressed in terms of the RHS of 
(3) as 


d e 

fF = Ac + Bu®™ + 7 aja) + S% BuO), &=0,1,2,...,N—1, 
j=l il=1 

(22) 


k-v1) k—-v2) k—va) 


where #*) = {x va! yee ol } and 
GM) = {ylk-wr) ylk-w2) lke) 1. Adapting (21) and (22) to (3) of the 


OCP yields the recurrence relation below: 


d e 
Bar) — FD 4 wyl®) < DS ajeh-es) — 5S > Byu®-’) (28) 
I=1 


iat 


where 6 = (I, + Ad) € R"*", w = Bd € R"™™, I, € R"™” (identity) and 
k =0,1,...,N-—1. Slotting the various values of & into (23) forms the linear 
inequality below: 

Az + Ba < Ex" + Fa" =C, (24) 
where A is a multi-diagonal matrix with nN rows and nN columns (ie., 
A € R'*%): B is also a diagonal matrix with nN rows and m(N + 1) 
columns (i.e., B € R'Y*x™N+)), and & is the concatenated vector of the 
state variable with dimension nV x 1 while u is the concatenated row vector 
of the control variable with dimension m(N + 1) x 1. The dimensions of 
the matrices EF, £”, F and a” are nN x n(va +1), n(vat+1) x 1, nN x mwe 
and mw- x 1, respectively. The respective structures of the concatenated 
vector-matrices are represented below: 


—I, 0 (Zee. Setar Nes atee. aeake ands 0 a) 
6 -I, 0 O : 
git) 
0 6 -I, 0 O 
p(v2—e1) 
ze ay 0 @ —-I, 0 : 
gp (va-Va-1) 
Qd ea QA eee 0 0 —Iyn (vet) 
0 da ay: O 6 -In 
Bete Lo. 0 (N-1) 
0 O Qa: a 0 ¢ —-I, (N) 
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w 0 0 --0 y) 
0w 0 0 
g(vitt) 
; Bie» Ow O y(e2—v1) 
Bu= ; 
Be-++ Br +++ Ow 0 
0 Be ; By ‘ y wet 1) 
0 0 Be By Ow QO u) 
—0 -a1 —as —Qd-1 —Qa (0) 
ate 0 gl) 
—Q\ —-As —Ag-1 —Qd 0 : 
meee glv2—vs) 
—-As —Agq-1 —Qd 0 
Ez = =“ 
—Qq-1 —Qd 0 
g(va—Va-1) 
—Ad 0 F 
Bette w(l—va) 
_ 0 g(-va) 
and 
— By — Bo —Bs —Be-1 —Be (0) 
Su 0 uD 
— By —Bs Os —Be-1 —Be 0 : 
5 yw2—ws) 
_ —B, aa —Be-1 —Be 0 
Fa" o= a ; 
—Be-1 —Be 0 
- ylwe—We-1) 
Be a3 0 
0 = yr we) 
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where the delay coefficients a and 3; are defined below as 


_ _ —dap, 1 < k < d, 

[ax] = { 0 elsewhere (25) 
and 

2.) — —OBr, 1< k < e€, 

[Px] = { 0 elsewhere. (26) 


However, the entries [@;;] € A, [bij] € B, [@,;] € E, and [fij] € F of the 
various matrix structures are described below as follows: 


In, t= 5 1<i<QN, 

_, J 6, jai-l 2<i<N, 

a ee ee ere ee sae ee 
0 elsewhere, 


(27) 


7 WwW, t=j LSi<QN, 
lbl=4< Pea j=i—w, I+ mpi nN, &=1,2,...(6, (28) 
0 elsewhere, 


-6, i=j=1, 

[é;;] = —AQp,7 =UR+1—-1 Lice, b= 1,2, sc0yd, (29) 
0 elsewhere, 

= Oy t=j=l1, 

[fiul= 4 Be j=wet1l—i 1<i<w,, k=1,2,...,e, (30) 


0 elsewhere. 


It is imperative to note that the real symmetric and positive definite prop- 
erties of the discretized matrices P,Q,A, and B affect the well-posedness 
of the algorithm and the selection of the optimal (relaxation and penalty) 
parameters; hence the bases for the assumption. 


3.4 Splitting method 


The D-R splitting technique, also known as ADMM, is a special case of the 
splitting technique and can be traced back to work done on monotone opera- 
tors and operator splitting methods into the general Hilbert spaces. Eckstein 
and Ferris [7] showed in turn that DR splitting is also a special case of the 
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proximal point algorithm applicable to optimal control systems. Suppose 
there exist two separable closed convex functions (or proximal operators) 
f(x) and g(a) that are non-expansive and not necessarily symmetric. Then 
the iterates of the D-R splitting algorithm for min f(x) + g(x) starting with 
any C° and repeating for k = 0,1,... are stated as 


iar ee Se + prox, (2«* =(\—2* 
ol <— prox, (4), 
where «* converges to the solution of 0 € Of(x) + Og(x), if the solution 
exists. The D-R Algorithm can be re-structured by introducing a new variable 
w* = x* — ¢* if when substituted into the D-R splitting Algorithm above, it 
yields its equivalent iterates expressed below as 


ae Prox, (a* + w*), 
ee oe prox -(u*t" —w*), 
wet wh 4 gktl — yb 


where x = prox;(u—w) satisfies 0 € Of (x) + Og(x). However, the equivalent 
D-R Algorithm can be amended by introducing the Nesterov-type accelerated 
variants or relaxation factor a € [0,2] in order to improve the rate of con- 
vergence of the iterates. The relaxed iterates is arrived at by replacing ut! 
with au*+! + (1—a)zx* in the sense of Nesterov [16]. The derived algorithm 
is stated as 


Initialize at any ¢°, x°, fixed a and repeat for k =0,1,..., 
ae prox, (* +w*), 
gh tt ¢ prox -(au*** +(1—a)z* — w*), 

k+1 ky gkt+l_ (y 


where it is an over-relaxation for 1 < a < 2, under-relaxation for 0 < 
a < 1 and reduces to its non-relaxed form if a = 1. The extension of 
the splitting method in the computations of the updates of the dual and 
adjoint variables of a separable convex constrained problem of the form 
min f(2%1) + g(x) s.t. h(a, 22) = 0 requires a sequential minimization of 
the augmented Lagrangian functional L,(x1, x2, ¢); hence the name alternat- 
ing direction method of multipliers (ADMM). 

The three steps iterates (updates) for the ADMM algorithm are as follows: 


ah? <— argmin (f(e1) + ¢Th(ar, 28) + £ || h(ar,2$) 13), (81) 
Ly 
xh*! <— arg min (g(a) +67 ACh, 92) + £ || awh, 22) 3), (82) 


Ch <— CB + ph(att*, 2577), (33) 
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where p > 0 is the penalty parameter and oe = axkt!_(1-a)ak for any ac- 
celerator variant a € [0,2] puts the algorithm in a relaxed mode. The ADMM 
is a more general method for handling all kinds of convex optimization prob- 
lems which includes the ¢;— regularization of the for min }°, fi(w) + ||w|, 
which other traditional methods such as gradient descent, Newton, an so 
on cannot handle. The ease with which the ADMM can be implemented 
in parallel for large-scaled data at a faster rate gives it an edge over other 
algorithms. 


3.5 Analysis on bounded control problems 


In a general box constrained nonlinear programming problem 
min{f(z)|¢EeQCR"}, (34) 


where 2 = {x € R”|l < & < uw} is a feasible space and a closed convex 
subset of R”, f(a) and a strongly convex quadratic function that is twice 
continuously differentiable in R”, J and u are the lower and upper bounds, 
respectively. Carreira-Perpindan [5] and Combettes and Pesquet [6] extended 
the theory of proximal operator to a proximal bound-constrained quadratic 
program (QP) (34) with the Lagrange formulation presented as 


Lp(a, z,€) = f(x) + g(z) + le—z2+é|2 st. x =z, (35) 


where z € R” is the slack (auxiliary) variable, € = 4 € R” is the scaled dual 
(adjoint) vector, p > 0 is the penalty parameter, and g(z) = § || z—v ||3 is 
the convex term (or indicator function for the nonnegative orthant). 

The proximal M-ADMM formulation is then expressed in the form 


a*+1 <— prox, p(2* — ¢#) = argmin (f(x) + £ || 2— z* +e" 3), (36) 
xz 
ae ; . (L p 
z+? 4 — prox, o(a*+1 — g) = argmin (£ || z—v 3) + & | att — 2 +e IB) 
st. l<z<u, (37) 


ght go ek 4 aktt -_ gktl. (38) 


However, there exist three possible cases of the location of the parabola 
vertex of the objective functional upon which the optimal solution is deter- 
mined. See sketches below for illustration. Each of the cases indicated in 
Figures 1, 2, and 3 represents the locations of the parabola vertex z7 within 
the bounds, below the lower bound, and above the upper bound, respectively, 
while Figure 4 represents the median of the three points. Therefore, the op- 
timal solution defined on the feasible space is the median of each of the ith 
component (7 = 1,2,...,N) of the three vectors; lower bound J, upper bound 
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LY wa 
a u. 7 U. 
li ae i 4j 1, i 


Figure 1: Vertex within bounds: 1; < Figure 2: Vertex below lower bound: 2; < 
yi Lu iF 


M( lis ZU; ) 


Figure 3: Vertex above upper bound: Figure 4: Median points: M(li, zi, wi) 
Ui S 2% 


u, and the parabola vertex z*. The choice of the next iterate z is determined 
as follows: 


M (27, li, ui) = Min |[Max(z}, l:), us| = Min|{l;, ui = li, af< lk, 

A= M(,2*,u) = Min|Max(I;, 27), us| — Min|[z¥, us] =e, 4 <2 < uy, 
M (lj, ui, 27) = Min|Max(l;, ui), 2] — Min|[u;, 27] =U, 2 > Uj. 

(39) 


3.6 M-ADMM algorithm formulation 


The original ADMM formulation is modified to accommodate both the state 
and control variables before imposing the Karush-Khun—Tucker (KKT) op- 
timality conditions. The Gauss-seidel accelerator variant is introduced to 
improve the convergence properties of the algorithm as illustrated in the 
convergence analysis below. The re-formulated compact convex quadratic 
optimization problem is then stated as 
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1 1 
MinJ(, au) = ; ET PE + xu Qt +R (40) 
s.t. G(x,u) = A+ BUu-C <0, (41) 
VSU<G, (42) 
ra > 0, (43) 
“vy+ua> 0, (44) 
where 
B= (ec a ™%) ER, Ga (wu, 2. ul) eR™ATY 


A E Rea, C E ine. B E Re sels) 
E Rey xn and a) E Resa 


while ¥ = (1/VN + 17, Ede ND yD), wg (+t) e RAT) and 
a = (1/VN + 1)(0? 0, oP oo N+) © R™N+D for all 
7, € y and o; € o, respectively. In addition, P and Q are the real, symmetric, 
and positive-definite matrix operators. 


3.7 Proximal-bound formulation 


In formulating the convez-proximal bound constrained program for the dis- 
cretized multi-delay OCP in (40)-(42), the Carreira-Perpind approach [5] 
is applied. Suppose that the control bound (42) is re-formulated into two 
equations, ¢ —-u > 0 and —y+%> 0 such that 


m 

0<5 | (o — 7) — (v1 + v2) [13 
Says eS 

y v 

<5 \|e-a-m [3+ 


Fl=y+@-@ 18 +5 Ily-0ll3, (45) 
where g(y) = § || y—v ||3 is a control bound function defined for u > y 
and 7 < y < a; while v is the slack variable replacing the inequalities and yu 
is the penalty parameter. Therefore, the associated augmented Lagrangian 
functional for the quadratic formulation is expressed as 


ae Loe Hh ae LL 
Lipy,p2(, 4, ZV, Ys E1, €2, b) =58) PEt su Qu + R+ 2 | y—v II3 


O || Az+ Ba C+z24+& |[2 
eae 
2 


|a-—yt & |I5 
st. z,v>Oand7<y<a, (46) 
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with the M-ADMM iterations stated as 


. 1 i 
gett g_ argmin {527 Pa = || Az + Buk —C +4 2* + € (3 1, 
xz 


1 ' 
att argmin { 5u7Qu + = || Az’** 4 Ba-C42* |2 4 2 | a—y* + €% [2 \, 
uu 


ghtl y argmin { © || Az*+1 + Bart! C44 ef |[2 \ s.t. z >0, 
z 


petty argmin { lly—v ||3 \ s.t.v>0, 
Vv 


+ A Le p22, _ _ 
yh? <— argmin {£ ly—v 2 +2 a! —y +e Bb st. 7s usa, 


eptl = ef + (Ag*t? + Ba*t? c+ 2**), 


eft = eb + (ah — yh), 


where & and tare the primal variables, \ is the Lagrange multiplier, p1, p2 > 
0 are the penalty parameters, ||-||2 is the euclidean (spectral) norm of a vector 
(matrix) argument, €; = A;/p:i, 1 = 1,2 are the scaled dual vectors, y is the 
introduced auxiliary (slack) vector, and /,(z) is the indicator function for the 
non-negative orthants defined as 1+(y) = § || y — v || and oo otherwise. 


3.8 Derivations of the update formulas 


Applying the KKT optimality conditions on the augmented Lagrangian func- 
tional (46) for the sequential minimization of the variables requires updating 
all the critical variables as indicated in the ADMM iterations. Moreover 
the update formulas of the adjoint variables will be accelerated using the 
Gauss-Seidel relaxation factor a € [0,2] as clearly illustrated in the works 
of Nesterov [16] and Ghadimi et al. [8]; though this work considers an over- 
relaxation where a € (1, 2]. 


State-update - Z: Setting 0;L(Z,-) = 0 yields 
g*t1 — _9,(P +p, A? A)1AT(Ba* —C+2* +f) (Z-— update), (47) 
where (P + pA’ A)~' is invertible since P is real, symmetric, and positive 


definite and A’ A is positive semi-definite sequel to Definitions 1, 2, and 3. 


Control-update - w : Setting 0, L(z**!, a, -*) = 0 and replacing Az**! with 
h(z**1, ak, a) = wAz**1 — (1—a)(Bu* —C + 2*) to relax the algorithm with 
the relaxation factor yield 


a*t! = —p(Q + pB" B)'B? [a(Az*t! —C + z*) — (1-0) Bu* + v*]. (48) 
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Auxiliary (Slack)-update - Z : Setting 0,L(-**1, z,-*) = 0 and relaxing 
the algorithm by replacing Az*+! with h(z*t!, a*, a) yields 


2* =—a(Az*t! + Bak! — C)— (1-a)[B(a* — a*) — 2*) — eF. (49) 
Therefore, 


gktl =max {0, zt (50) 


=max {0,-a(Az*+! + Bak+ _ CO) — (1-0) [B(akt! — a) — 2] - et}. 


Updating the parabola vertex - y : Setting 0,L(-"t',y,-*) = 0 in the 
derivation of the parabola vertex y in the box constraint yields 


OyL(***,y,-") = uly — v) — po(i—y + &2) =0 (51) 
k+1 oktl 1 ¢k 
y* = MU + p2(u + &) s.t. 7 < y* 2 a, (52) 
(+ pa) 
where for each component y* of y* = (y1,y2,---,YNn), the minimum of the 


upward parabola is defined in the sub-interval [4;, 7] such that the solution is 
the median [4j;, y7, o;] of 7,4; and each component y* of the parabola vertex 
derived above is specifically defined as 


. Beet! + po(azt + &&) 


y;, = s.t. Fi < yf <o; and element-wise for all 7. 
(u + p2) 
(53) 
Therefore, 
k+1 —k+1 k 
yh tt = Min y, LU + p2(t +8) 5 ; (54) 


(p2 + 11) 
= Min[Max(7, y*), a], 


with the steps involving element-wise thresholding operations. Considering 
the extension of the above concept of proximal bound operator on the OCP, 
let the state variable 7 € R% belong to a set of admissible state functions 
X (written Z € X C R©) and let the control variable a € R(Y+) belong 
to the set of admissible control functions U (written @ © U € R™). We 
then assume that v7 = (#*,7) and v\%) = (z*,a) are the points on the 
state-control coordinates (%,%) for which the control is bounded below and 
above by 7 € Rt) and & € Rt, respectively. If there exists any 
point on the ith component of the control vector (¢ = 1,2,...(N +1)), then 
yi € [7i, Fi] is such that yi) = (Z%, y;) is the point on the (2, i) trajectory 
that prescribes the parabola vertex of the objective function f(Z, i) at its 
minimum. Therefore, the median M of v‘%), v(%) and the parabola vertex 
v4) are defined as follows: 
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Min [Max(v%), v@)), va] = yi), yi) € [v%), ya], 
(at a) = Min [Max(v@), o(%)), ya] = y%), yi) < yl), 
Min [Max(v), yi), yo] = pF), yi) > yi), 
(55) 
where % = £* and u = y; are the optimal state and control variables, respec- 


tively, within the feasible space. Therefore, at Z = %‘*+!), the G-update is 
defined as 


(z+), yt) =M [vO v2), v9] , element-wise for all i € {1,2,...,(N4+1)}, 
(56) 
where v\¥) = (@*,y*) = (@%, yt") such that 
k+1 _ port? ob pa(urt* ae &k,) <* _ —(k+1) 
YU; 3 x; x; ; 
(p2 + H) 


ui; =yi, for alli € {1,2,...,(N+1)}. 

(57) 
For further illustration on the analysis of the proximal bound and its effects 
on the objective functional value, we then consider the figures below. 
Suppose that (Z*, y*) is the vertex parabola for which f(%, a) is at minimum. 
Then 


F(z, y*) < F(Z, @), st. 2 © X, forallieUu and ye 7, a]. 


Assume that there exists a bound 7, a] C U on the control variable, sim- 
ply called bounded (box) constraint, for which the objective functional value 
F(,t) is at minimum; then in the process of the algorithm, the three ex- 
pected cases are stated as 


F(z*,y*), for all a € [7, |] CU Case I 
F(z ,y*)=4¢ F(#,7), forali<7CU Case II (58) 
F(z#*,o), forali>aCuU Case III 


CaseI: y<u<oa 

In Figure 5, the vertex point (%*,&%”) that falls within the control bound is 
same as the optimal state-control trajectories (z*,u*) such that F(Z, %*) = 
F(z,u"). The functional Q within the region of feasibility is then defined as 


Q* = {F(z*,u”)|z* © X and ue |[4, a]. (59) 


Case II: u<7 

In Figure 6, the vertex point (*, a") that falls below the lower control bound 
is same as the optimal state-control trajectories (<*,7) such that F(Z, u*) = 
F(z,7). The functional Q within the region of feasibility is then defined as 


O* = {F(z*,7)|z* © X andu<¥_ forall ae U} (60) 
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fw) 


s\ 


control— 


Figure 5: Control within bounds: y < Figure 6: control below lower bound: 
ue <o a <7 


Case III: uw>a 

In Figure 7, the vertex point (%*,t”) falls above the upper control bound 
with optimal state-control trajectories given as (Z*,@) such that F'(%,i*) = 
F'(#,a). The functional 2 within the region of feasibility is then defined as 


O* = {F(z*,7) |e eX andti>a  forallue UV}. (61) 


Updating the scaled-dual variables - €,,  : It yields 
RH ek 4 (Apet1 4 Baktl — C4 2k), (62) 


and 
i ee i ie (63) 
3.9 Primal-dual convergence 


The convergence of the M-ADMM for the bounded control in terms of its 
primal-dual feasibility requires deriving the primal-dual residuals and their 
behaviors for large iterations. 
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Figure 7: Control above upper bound: 
a<@ 


Theorem 3. Given the bounded control problem (40)—(42) with the mul- 
tipliers for the linear and bounded control constraints given as A; and Az, 
respectively, then there exists a dual residual s*t1 = p,BT[z*t! — z*] — 
p2 aoe - y*| that converges to zero for given penalty parameters p; and p2 
for the linear and control bounds, respectively. 


Proof. Given the objective function p* = f(a*) + g(u*), linear inequality 
constraint Ax + Bu < C and y = o—7, the associated Lagrangian with slack 
z is stated as 


L(, U, 2, A1, A2 : pi, p2) = f(%) + g(a) Nj (Az + Bu C +z) 
+5 || Ab + Ba-C +213 +5 lly-o lk 


- P2 = 
M(u—y) + F | a—y IB 


Applying the optimality conditions KKT to the Lagrangian function above 
yields 
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Ag(u) + BT [Ak + p\(Az*t! + Bu— C+ 2*)] + AF + po(a—y*) =0, 
Og(ukt) + BT [AT + pi (AB*t + Bakt? —C 4 2**1)] +28 
ee 


k+1 
V4 


+ po (G*** — y*¥*) —p, BP ZA) + py BE z* + pay**? — poy* = 0, 
a, 


k+1 
"3 


Og(ukt!) + BTM + pirf tt |] +A$ + port? 
{jV——-_“_‘ n_—>¥/-@"/[—“— 


AP ABH 
= Br [geet _ 2" — ps gt = a | 
—E EES Ss 
gktl 


where the primal residuals are expressed as 


htt = Agtt) + Bukt) — C+ 2**1, (64) 
pee = gett _ yer, (65) 


If «*+1 = x* minimizes the convex function, then 


Oq(a*) + BT +5 =0 (66) 


and setting u*+! = y*+1 


as 


completes the proof by expressing the dual residual 


gkt1 = p, BT ae _ 2" — po [eet _ y*] . (67) 


Let the state-control variables be concatenated as w = (Z, u) with dimen- 
sion nN +m(N +1) x1, A= [A,B] with dimension nN x (nN +mN +m), 


and P = wivaerel with dimension (nN +mN +m) x (nN +mN +m). 


This yields a re-formulated augmented problem of the form 
wo PH+R st. AW<C, (68) 


which results into Theorem 4, as an extension of Ghadimi et al. [8]. 
Theorem 4. Given a convex quadratic programming (QP) problem r 
min 50? PH+q' w such that Aw < 6, where P € R"™", WER", qER",A 
R™*”, and b € R™, then the optimal step-size for the QP is 


-1 


i= / Amin (AP-1AP) max (P-L) (69) 


and the convergence factor is 
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(70) 


For a € (1,2],0 < ta <1, P is symmetric and positive definite and 
Xmin(AP~!A?) and Amax(AP~'A™) are the minimum and maximum eigen- 
values of the matrix AP~!A’, respectively. 


3.10 The Error dynamics of the relaxed Algorithm 


Let us consider the augmented convex OCP in (68) with ¢2-regularized esti- 
mation, with the regularization parameters 6 > 0, as stated below: 


P 6 r 
a? Pa + R+ 5 llalla st. Adb+2-C=0. (71) 


The resulting relaxed ADMM iterations take the form 


ort — — (P + pA A) = AT (\* + p(z* — C)) (72) 
. [>" +p (a(Aatt? —C)-(1- a)z*)| 73) 
d+ p 


vetl =F 4 5 (a(Awt? +2841 _ ©) 4 (1—a)(z#t! — z*)) _ (74) 


Considering the convergence of the relaxed iterations, (72)—(74), it is imper- 
ative to note that the relaxed iterations return the standard iterations of the 
convex control problem by setting a = 1 and 6 = 0 with the augmented vari- 
able w*+1 = [@*+! Gk+1|7_ In characterizing the convergence of the relaxed 
iterations, we then analyze the error dynamics of the ADMM to know how 
the errors associated with w* or z* vanish. Making \* subject of formula in 
(73) and inserting into (74) yield \*+! = —d5z**!. Inserting the w-update 
into the z-update and using the fact that A* = —dz*, we arrive at 


1 eed zk n 
k+l _ = = T A\—-1 {AT k 
z = 55, [+e (a a)I +a(p —5)A(P + pATA)-1A NE 
— 
E 
ap Was AT A\—-1 AT 
=P | pAtP + pAT Ay-1 AP — 1] ©. 
5p PAP + 047 A) G (75) 
For z = 2* a fixed-point of (75), that is, 
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*_ myx _ OP as AT A\-1 AT _ 
z= Ez 549 PAP + 0A A) A’ -1/C, 


then the dual error e*+! = z*+1 — z* evolves as e&+! = Ee* and z*t1 — z* = 


E(z* — z*) such that the z— update of the relaxed iterations above converge 
if and only if the spectral radius of the error matrix EF in the above linear 
iterations is less than one. The eigenvalues of F can then be written as 


ap(A;(P) + 6) 


RAO, ri P= A ’ 
Se EG CAPD) 


(76) 


for p, 6 and A,(P) € Ry (ie., positive) provided a (1, oer P ters) ) 
To ascertain the upper bounds p and a that guarantee an improvement with 
strictly smaller convergence factor compared to that of the classical ADMM 


iterations, we take the derivatives a9 with respect to the parameters p and 


a given by (p*,a*) =argmin max |éa(a. p,A:(P))| while that for the lower 


bound of €% is & = max lea(p", a*,\;(P))|. This yields the jointly optimal 
step-size and gonvensenod factor expressed as (p* oh *), where the detail proof 
of the convergence analysis can be seen in [8]. We hea conclude that the over- 


relaxed iterations are guaranteed to converge faster for all a € (1,2]. The 
developed M-ADMM algorithm is stated below; see Appendix for flowchart. 


3.11 Algorithm: M-ADMM for proximal 
bound-constrained program 


Input parameters and operators p, a, A, B, C, PM, edual 
Initialize 7°, w°, 2°, £2 = A9/o1, €8 =AP/p2 

Set P>0O, Q>0 (symmetric and positive definite) 

For k=0,1,2,...,do 


Compute ght using (47) 

Compute ue ee (48) 

Compute 1 using (50) 

Compute ee +1 using (54) and set af*t = y**t,  element-wise 
for all 2 


Update éft! and é$t! using (62) and (63), respectively. 
Compute — ||r**1||3 and ||s**+||? using (64) and (67), respectively. 
if |jre+2)/2 < ePTiM and \|s*+2)|2 - -dual 

| stop and return (output) 

else 


| repeat process 
end (If) 
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end (for) 


Return 7*t1, gktl, gktl ye, gktl gk+l 


’ u > SL » S2 


3.12 Computational Complexity 


Iteration complexity: Each step of the M-ADMM updates (iterations) is O(V) 
with large computational efforts required. However, catching the Cholesky 
factorization of the matrices (P + p;A7 A) and (Q + p1B?B + p2I) in the 
form RR? (with R the upper triangular) will make the computation more 
effective, provided P and Q are real symmetric and positive definite. If P or 
Q is dense, then the computation of the Cholesky factor is a one-time cost of 
O (3.N*) while the linear system is O(N”) by solving the two triangular sys- 
tems. However, if P or Q is large and (sufficiently) sparse, then computing 
the Cholesky factor and solving the linear system are both possibly O(N). 
The solution can be obtained in few iterations of the M-ADMM algorithm 
when the CGM or L-BFGS iterative solver is deployed with a warm-start to 
carry out faster convergence of the updates. In practice, the warm-start in 
which the values of z = v and €; = 2 = 0 are usually adopted. 

Termination criteria: The earlier derived stopping criteria, in (64)—(67), can 
be deployed with a tolerance range of 10~* (for k = 3,...,6). Moreover, 
the speed at which the M-ADMM converges depends significantly on the 
quadratic penalty parameters by Boyd et al. in [3]. However, the primal 
ePrim > 0 and dual «?"*! > 0 tolerances are selected as the termination 
(stopping) criteria for the convergence of the M-ADMM, and their choices 
depend on the relative «"® and absolute €%* criteria on account that the 02 
norms are in R” and R™, respectively. The selected primal and dual toler- 
ances are so small such that the algorithm converges (terminates) whenever 
I|rFt4]]o < €Pr'™ and ||s**1\|2 < e?4. Usually in literature, the values, 
eT! — 1073 and €%* = 107%, are used as reasonable stopping criteria for 
the ADMM algorithms and can be reduced to improve the accuracy. Stated 
below are basic computations in literature as in [3]. 


ePrim = s/ne* + max{||Ax***||p, ||Bu***||2, |[2***]]2,||Cll2} (77) 
cDual = a/mers 4 ce" |loullo (78) 
4 Numerical experiments 


In this section, the algorithm was implemented with numerical examples on 
a Core i3 CPU. The quasi-Newton solver is interfaced with the M-ADMM 
in obtaining the optimizer. The convergence factor of ADMM for varying 
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tolerances and penalty parameters of the linear constraint is later illustrated. 
These examples demonstrate that the M-ADMM method converges for an 
over-relaxed factor a* € [1.8, 2.0] for faster convergence. 


Example 1. Consider the one-dimensional convex multi-delay bounded 
OCP: 


Miceniins Pe a = : | [202(t) + u2(t)]at 


subject to a(t) < —2e(t) + 3u(t) + 2x(t — 0.1) + x(t — 0.2) 
—a(t — 0.5) + u(t — 0.1) — 2u(t — 0.3) 
x(0) = 1, t € [0,1], 
a(t) = 2é7 +41, t € [-0.5, 0], 
u(t) = 3t+2, t € [—0.3, 0], 
0< u(t) < 2. 


The discretized positive definite matrices P(> 0) and Q(> 0) for the objective 
functional, using 6 = 0.1 for purpose of illustration, are, respectively, given 
as 


0.0889 0 --- «=. 0 0.0444 0 -- «0 
0 0.04440 -. : 0 0.0222 0 
P= and Q = 
"0.0889 0 +. 9.9222 0 
O ses ees 0 0.0222 0 ike oes Q) Omid 


While the constraint coefficients are as given in subsection 3.3 using the 

optimal relaxation parameter a* = 2.0. Table 1 demonstrates the various 
effects of changing relaxation parameters (factors) on the convergence and 
results of the M-ADMM algorithm measured by the iterative cycles. The 
optimal relaxation factor is observed to be a* = 1.90. 
Figures 8 and 9 show the effects of the varying penalty parameters pz on the 
rate of convergence of the M-ADMM algorithm for increasing relaxation fac- 
tors a. Clearly, the lines overlap as the algorithm progresses. The optimum 
values were also obtained at a* = 1.90, as shown in the simulation, indicating 
an increase in the rate of convergence. Table 2 shows the results when the 
algorithm was ran for various values of the step-length 6 and tolerances for 
fixed values of p1, p2, 93, and p. 

Figure 10 illustrates primal and dual convergences of the M-ADMM al- 
gorithm using the optimal over-relaxation factor a = 2.00 while Figure 11 
demonstrates the responses of the state-control variables and the objective 
values for the chosen tolerance (Tol.=10~3). The objective values J* = 0.3824 
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Tolerances 1073 10-4 10-5 10-° 
No. of iterations (k) (k) (k) (k) 
a = 1.00 66 128 193 259 
a = 1.20 55 107 160 216 
a = 1.40 A7 98 138 189 
a = 1.60 42 80 122 170 
a = 1.80 37 71 109 155 
a* = 2.00 34 69 130 204 
Table 1: Effects of varying relaxation factors on fixed penalty parameters (ju = 0.3, 91 = 
0.1, p2 = 0.2, p3 = 0.3:5 =0.1) 
Step-lengths 6 = 0.10 6 = 0.05 6 = 0.01 
Tolerances | 10-7 107% i0-? -10-* 10-340? 
||| 2.3788 2.3919 2.8830 2.9223 3.2970 3.6359 
||] | 0.8937 0.8933 0.9271 0.9263 0.9812 0.9796 
|| J * || 0.3827 0.3900 0.3968 0.3982 0.3623 0.3730 
oh 0.5827 0.5827 0.0528 0.0528 0.0471 0.0471 
PCAL 0.0557 0.0557 0.7379 0.7379 0.9383 0.9383 
k 16 14 25 34 69 
Time (secs) 0.0504 0.05623 0.0886 00305 0.3778 0.2057 
Table 2: Summary of results at optimum relaxation parameter (u = 0.3, p1 = 0.1, p2 


0.2, p3 = 0.3, a* = 2.0) 
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Figure 8: Effects of varying a, pon conv. Figure 9: Effects of varying a, Tol. on 
as ||r®||2 < e¢¥@! and ||s*||2 < primal conv. as ||r*||2 < «4&4! and ||s*||2 < 
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0.2 0.5 ————— 
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number of iterations number of iterations 


Figure 10: Primal-Dual Convergence: Figure 11: Optimal trajectories: State, 
|r*||2 Sa a \|d*||3 >0 control, and objective 


were arrived at with the respective optimal state and control values 2.3788 
and 0.8937 € [0, 2], respectively. 


Example 2. Considering the generalized multi-delay OCP with bounded 
control: 


0.5 
MinJ (a, w) =| [et + e102 +05 + uy + uru2 + u3|dt (79) 
0 


s.t. 21 (t) <2a1(t) + wo(t) + wi (t) + 3ue(t) + v1 (t — 0.1) + xo(t — 0.1) 
+ 2x1(t — 0.2) + ro(t — 0.2) — a1 (t — 0.3) + 2ui(t — 0.1) 
+ u9(t — 0.1) + ui (t — 0.3) + 2ue(t — 0.2), 
Lo(t) <a, (t) — u(t) + 2ua(t) — x(t — 0.1) + x(t — 0.2) 
+ 2%o(t — 0.2) — a1 (t — 0.3) + wa(t — 0.1) + ur (t — 0.2) 


+ 3u2(t — 0.2), t € [0,0.5] (80) 
«(0) =(1, 1), (81) 
a(t) =(2t + 1,4? +1), —0.3<t< 0, (82) 
u(t) =(2,2 +t), =(2 272 6, (83) 
(1,1) < u(t) < (2,3), (84) 
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where x(t) € R? and u(t) € R?. 
The problem is re-formulated into the format in (1)—(7) to make it amenable 
to the M-ADMM algorithm. The coefficient matrices are defined below. 


p= [fi]. = [28]. 4= [2]. 2=[43]-2= [=f 


a= |i g)e2=([7a]e2= [270]. = [or] ane ee= [73] 


The coefficient matrices, P and Q, of the quadratic functional are sym- 
metric, invertible (non-singular), and positive definite (ie., P > 0 and 
Q > 0) since their respective eigenvalues (Ay = 1.585, Ag = 4.4142) and 
(A, = 1, A2 = 3) are all positive. This ensures that the operators are well- 
posed for the algorithm. The choice of the step-length 6 = 0.1 was made 
with N = 5 such that the coefficient matrices P and Q of the state and 
control variables for the continuous-time problem yield the block matrices 
P € R'*!0 and Q € R!*!?, respectively. The structures of their block 
matrix operators, are as defined in (3.2). However, the recurrence relation 
from the discretization of the constraint, as derived in (23), yields 


3 2 
Oar) — FAD 4 wry (®) = 6° avja'h-ea) — 5S Brum), (85) 
j=l i=1 


where v = 13/0 = 3, w = q2/d = 2 and the various associated matrices are 
defined below: 


O=1,4A5= bee aon pee | 0.1000 al | 


0.1000 1.0000 —0.1000 0.2000 


with the formulated block-matrices A, B, and C are as defined in subsection 
3.3. Table 3 shows the results in the first phase of the algorithm. The 
M-ADMM was ran for various values of the relaxation parameter a evenly 
spaced between 1.0 and 2.0 with a step-size of 0.05 for a fixed value of the 
penalty parameters (p; = 0.1, po = 0.2, ps = 0.3, a* = 2.0) and tolerance 
(Tol = 10~%) to ascertain the optimum relaxation factor (a*) using the 
number of iterations as the measure of convergence. 

Table 4 also demonstrates the effects of some of the selected relaxation pa- 
rameters (factors) on the optimum values of the optimization problem as the 
tolerances decrease geometrically for fixed values of the penalty parameters. 
Figures 12 and 13 show the rate of convergence of M-ADMM Algorithm using 
the primal and dual residuals for increasing values of the iteration. Clearly, 
both lines converge to zero as the algorithm progresses for the varying values 
of a and p2 as shown in Figure 14, where the average elapsed time per/cycle 
falls within a second when run on an Intel computer (Inspiron 15 Core 13). 
The optimum values were obtained at a* = 2.0, as shown in the number of 
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Step-lengths 6 = 0.10 6 = 0.05 6 =0.01 

Tolerances 10-3 10-4 10-3 10-7 10-3 10-* 
Tal] 14.5367 14.1757 0.138099 0.135782 0.138099 0.135782 

|||] 2.7139 2.72038 0.15 144 0.22 164 

|| 7 * || 18.3905 18.6315 0.16 140 0.22 164 

& 0.6910 0.6910 0.16 140 0.22 164 

PCAL 0.1612 0.1612 0.16 140 0.22 164 

k 117 138 0.16 140 0.22 164 

Time (secs) 0.2678 0.3852 0.16 140 0.22 164 


Table 3: Summary of results at optimum relaxation parameter 
(p1 = 0.1, po = 0.2, p3 = 0.3, a* = 2.0) 


Tolerances 17? 10-4 10% 10-4 
No. of iterations (k) (k) (k) (k) 
a= 1.00 239 467 703 940 

a = 1.20 199 389 585 786 

a= 1.40 175 345 523 700 

a = 1.60 155 307 464 622 

a = 1.80 140 277 419 562 

a* = 2.00 127 253 383 514 


Table 4: Effects of varying relaxation factors on fixed penalty parameters (p1 = 
0.1, p2 = 0.2, pg =0.3:6=0.1) 


iterations per cycle, indicating a decrease in the rate of convergence. The 
M-ADMM ran for various values of the step-length 6 geometrically spaced 
between 0.010 and 0.10 with a constant factor of 0.5 using the tolerances of 


10-3 and 10-4. Figure 14 further illustrates the faster rate of convergence 
a5 25 
: convergence for —ilnil, (Primal) 
pal = 0.1, w= ||s||2 (dual) 2 —l[s||, (dual) 
g , p2 = 0.2 and g 
Zo.3 p3 = 0.3 at B15 
8 a* = 2.0. 3 | 
a, prmardualsing co Oe and 
wo p,=0.3 at the opt. relaxation 
£ parameter «.=2.0. 
a 
05 
% 0 100 150 0 
5 ; : 5 0 20 40 60 80 100 120 
number of iterations number of iterations 


Figure 12: Dual convergence: |{d*||3 +0, Figure 13:  Primal-dual convergence: 
pand a* = 2.0 lr*||2 > r*, ||d*||2 +0 


of the M-—ADMM with the optimal over-relaxation factor a = 2.00 compared 
to the non-relaxed ADMM for which a* = 1. However, Figure 15 demon- 
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Figure 14: Convergence for varying a, p2: Figure 15: Optimal trajectories: State, 
(p1 = 0.1, p3 = 0.3, a* = 2.0) control, and objective 


strates the responses of the state-control variables and the objective values to 
increasing iterations at chosen tolerance (Tol.=10~3), relaxation parameters 
(p, = 0.1, p2 = 0.2, p3 = 0.3) and optimal penalty parameter a* = 2.0. In 
addition, the objective values are observed to increase with increasing values 
of the state variables until the optimal objective value J* = 18.3905 was 
arrived at with the respective optimum state and control trajectories stated 
as 


» [41416] og ye = [1-000] rag [lle"ll] — [145367 
* = | 13.9343 “= | 2.5230 Ju*||| — | 2.7139 


and ||(1, 2)7|| < |Ju*|| < ||(2,3)7||. Applying the formulas in (69) and 
(70) yields the results for the calculated penalty parameter and conver- 
gence factor as p* = 0.1612 and zi; = 0.6910, respectively, where 
P € RON+2)x(4N42) and A = [A, B] € R2N*(4N42), 


5 Conclusion 


The ADMM for strictly convex QP with a particular focus on the QPs arising 
from discretized multi-delay OCP with bounded control was studied. A con- 
vergence proof and a method of analysis of the problem were proposed that 
allows the selection of the optimal relaxation parameter for the acceleration 
of the algorithm for varying penalty parameters of the augmented Lagrangian 
functional; while other parameters remain fixed. A major drawback of this 
approach is that, the penalty parameters for both the linear and bound con- 
straints cannot be derived simultaneously. In the future, we will extend the 
computational framework to include the effects of both varying parameters in 
the computational complexity of the algorithm. This research was limited to 
OCP with linear constraints of non-fractional order. However, in subsequent 
work, the scope of the work will be extended to convex constrained OCPs 
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with coupled and/or nonlinear differential constraints with both integer and 
(or) fractional order. 
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Appendix 


A. Derivation of P and Q 


By Simpson’s rule, if for each k = 1,2,..., x —1, f(t2x) appears in the 


term corresponding to the interval [t2,—2,t2x] and f(t2,-1) appears in the 
term corresponding to the interval [to,, tex42] for k = 1,2,..., x, then the 
approximation of the integral in (18), given that F(t, x(t), u(t)) = 2? Pa + 
u™ Qu, can be expressed as 


1 


T 
I= a F(t,a(t),u(t)) dt ~ 5 F(tp, 2), ul) (86) 
to k=1 
x71 = 
eee 2 (7 +2 57> fH 445 > pOk-v +7) 
~ 243 
k=1 k=1 
where 
FO = (we FP Px + (u™ )PQu, (87) 
fe") = (x 2h))P Py (2h) a (u2*) TP Qu l2h) (88) 
gee) = (g@e-)\t pyer-)) oe (ul2F#-D) Py (2k-V) and (89) 
FY) = (e&™)P Pe + (UM) TQuUM), (90) 


y-1 z 
ps ° (a) P (0) 4 E (a (24))P Py(2h) 4 Sa yPPalD 
k=1 k=1 
Nj 
6 176 26 % 
4 (g)yP Pal) de 5 | (ut Pu + = So (u@®)? Pu”) 
3 pa 3 
PS ons 5 
+ (uA) 7 Pus) + sy Pu), (91) 
k=1 


Expanding (91) for various values of & gives 
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46 
40D. fei ees 1 
- .\ fa 
0 2p- ry 
I wo (2)T Pa + (o,2,...,0) . 
i ali 
45D 9 
Qo ke aah ‘O $p oh’?) 
6 
2 Oe ies dake Bat site 0 
3@ AO) 
tC) i 1 
0 2Q uf) 
a (uo uf) _ ant ul) : * 2%Q é 
y=) 
45.Q Py uN?) 
OQ tus dw Ba 0 4Q 
expressed as 
_ ol — 1 = 
Iemin — 27 Pz+ —a7 Qu R, (92) 
zu 2 2 


where any elements p;; € P and Gj; € Q are described as 


2P, i =j(odd) 1<i<N-1, 
_ 22 =P, 4 = j(even) 2<i<N, 
Pul=) tp j2; 4=2N 

3 ? Sf); = ’ 

0 elsewhere 


and 


§Q, i=j, 1=1,N, 


a) = #2Q, i= j(even) 2245 N, 
@ 2Q, i= j(odd) 3<i<N, 
0 elsewhere. 


B. Flowchart for the M-ADMM Algorithm 


The diagram below (Figure 16) is the flowchart demonstrating the procedure 
of the M-ADMM algorithm in handling the discretized optimal control (or 
constrained optimization) problem with bounded (boxed) control. The left- 
side of the flowchart is the M-ADMM subroutine, which is the outer loop 
of Algorithm. However, the right-side of the flowchart is the optimization 
solver, which in this case is the ”quasi-Newton solver” of the low-memory 
BFGS type. 
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Figure 16: Flowchart for M-ADMM with bounded control 
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